MODULE PolynomialInterpolation
        USE precision
        IMPLICIT NONE
        PRIVATE
        PUBLIC :: Interpolate
CONTAINS
        !********************************************************************
        ! Given arrays xs(:),ys(:) of n x and y values, this SUBROUTINE will
        ! fit a degree n-1 polynomial P(x) to the data and return the value
        ! This SUBROUTINE is based on the Numerical Recipes SUBROUTINE polint
        ! Reference: Section 3.1 of Numerical Recipes
        !********************************************************************
        SUBROUTINE Interpolate(xs,ys,x,y,yerr)
                REAL(D),INTENT(IN)  :: xs(:),ys(:) ! array of x and y values
                REAL(D),INTENT(IN)  :: x
                REAL(D),INTENT(OUT) :: y,yerr
                INTEGER :: i,m,ns,n
                REAL(D)                     :: den,dif,dift,ho,hp,w
                REAL(D),DIMENSION(SIZE(xs)) :: C,bigD
                n=SIZE(xs)
                ns=1
                dif=abs(x-xs(1))
                closest: DO i=1,n ! Find index ns of closest table entry
                        dift=abs(x-xs(i))
                        if (dift.lt.dif) then
                                ns=i
                                dif=dift
                        endif
                        C(i)=ys(i) ! Initialize tableau of C's and D's
                        bigD(i)=ys(i)
                END DO closest
                y=ys(ns) ! Initial approximation to y
                ns=ns-1
                DO m=1,n-1         ! For each column of the tableau,
                        DO i=1,n-m ! loop over C's and D's and update them.
                                ho=xs(i)-x
                                hp=xs(i+m)-x
                                w=C(i+1)-bigD(i)
                                den=ho-hp
!***************************** This test will fail only if two input xs are
!***************************** identical within roundoff:
                                fail: IF(den.eq.0.0_D) THEN
                                        WRITE(0,*)'Failure in Interpolate'
                                        STOP
                                END IF fail
                                den=w/den
                                bigD(i)=hp*den
                                C(i)=ho*den
                        END DO
                        if (2*ns.lt.n-m)then
                                yerr=C(ns+1)
                        else
                                yerr=bigD(ns)
                                ns=ns-1
                        endif
                        y=y+yerr
                END DO
        END SUBROUTINE Interpolate
END MODULE PolynomialInterpolation
